#read ace files to estimate coverage for specific base

#see command_line_notes.txt for details on grouse1CO.txt 

readAceCO <- function(file.name){
	#read in line numbers from a grepped ace file
	#eg grep -n 'CO' file.ace > fileCO.txt
	AceCO <- read.table(file.name,stringsAsFactors=FALSE)
	names(AceCO) <- c("line","contig","length","no.reads","segments","orientation")
	AceCO$line <- sub(':CO','',AceCO$line)
	AceCO$line <- as.numeric(AceCO$line)
	AceCO
	}

#grouse1CO <- readAceCO('/Paterson/Datafiles/grouse/mapping/grouse1CO.txt')
#grouse2CO <- readAceCO('/Paterson/Datafiles/grouse/mapping/grouse2CO.txt')


readContigAce <- function(ace.file,AceCO,contig){
	#read lines in corresponding to a specific contig
	#from an ace file
	#AceCO object is returned by readAceCO function
	len.ace <- system('wc -l mapping/grouse1.ace',intern=T)
	len.ace <- sub('^[ ]*','',len.ace)
	len.ace <- strsplit(len.ace," ")[[1]][1]
	
	start.read <- AceCO$line[match(contig,AceCO$contig)]
	
	
	if(match(contig,AceCO$contig)==length(AceCO$line)){
		stop.read <- len.ace 
		}else{			
		stop.read <- AceCO$line[match(contig,AceCO$contig)+1]-1
		}
	tmp.cmd <- paste("sed -n \'",start.read,",",stop.read,"p\' ",ace.file,sep="")
	linesAce <- system(tmp.cmd,intern=T)
	linesAce
	}

#debug(readContigAce)
#tst <- readContigAce(ace.file='/Paterson/Datafiles/grouse/mapping/grouse1.ace',AceCO=grouse1CO,contig='contig01000')

#produce some complex list object
makeAceObject <- function(linesAce,read.info=FALSE){
	#creates a list object from a vector of character strings
	#generated by readContigsAce
	
	AceObject <- list()
	tmp.CO <- linesAce[grep('CO',linesAce)]
	AceObject[['CO']] <- strsplit(tmp.CO," ")[[1]][2:6]
	names(AceObject[['CO']]) <- c("contig","length","no.reads","segments","orientation")
	
	#padded consensus sequence
	AceObject[['padded.seq']] <- paste(linesAce[2:(grep("^BQ",linesAce)[1]-1)],collapse="")
	
	AceObject[['BQ']] <- as.numeric(strsplit(paste(linesAce[(grep('^BQ',linesAce)+1):(grep("^AF",linesAce)[1]-1)],collapse="")," ")[[1]])
	
	#reads and their position relative to padded consensus
	tmp.AF <- linesAce[grep('^AF',linesAce)]
	AceObject[['AF']] <- data.frame(
		read = sapply(X=tmp.AF,FUN=function(X){strsplit(X," ")[[1]][2]}),
		orientation = sapply(X=tmp.AF,FUN=function(X){strsplit(X," ")[[1]][3]}),
		start.pos = as.numeric(sapply(X=tmp.AF,FUN=function(X){strsplit(X," ")[[1]][4]})),
		stringsAsFactors=FALSE)
	
	#ignore BS
	
	#read lengths
	tmp.RD <- linesAce[grep('^RD',linesAce)]
	AceObject[['RD']] <- data.frame(
		read = sapply(X=tmp.RD,function(X)strsplit(X," ")[[1]][2]),
		length = as.numeric(sapply(X=tmp.RD,function(X)strsplit(X," ")[[1]][3])),
		stringsAsFactors=FALSE)
	
	if(read.info==TRUE){ #tried for mapCoverageQual
		tmp.read.seq <- list()
		for(i in 1:length(tmp.RD)){
			tmp.st <- grep('^RD',linesAce)[i]+1
			tmp.sp <- grep('^RD',linesAce)[i]+ceiling(AceObject[['RD']]$length[i]/50)
			tmp.read.seq[[i]] <- paste(linesAce[tmp.st:tmp.sp],collapse="")
			}
		names(tmp.read.seq)	<- substr(AceObject[['RD']]$read,1,14)
		}
	

		}
	AceObject[['read.seq']] <- tmp.read.seq
	#position of the read wrt padded consensus
	tmp.QA <- linesAce[grep('^QA',linesAce)]
	AceObject[['QA']] <- data.frame(
		read = AceObject[['RD']]$read,
		start.pos1 = as.numeric(sapply(X=tmp.QA,function(X)strsplit(X," ")[[1]][2])),
		stop.pos1 = as.numeric(sapply(X=tmp.QA,function(X)strsplit(X," ")[[1]][3])),
		start.pos.clear = as.numeric(sapply(X=tmp.QA,function(X)strsplit(X," ")[[1]][4])),
		stop.pos.clear = as.numeric(sapply(X=tmp.QA,function(X)strsplit(X," ")[[1]][5])),
		stringsAsFactors=FALSE)
		
	#watch out for 'contigs' in the reads (sub-contigs?)
	AceObject
	}
#tmp.ace <- makeAceObject(linesAce)

#structure of GS mapper ace files are unclear.
#Assume for now that one can ignore reads labelled as contigs
#within the RD slot

#calculate coverage for each base in a contig
mapCoverage <- function(AceObject){
	#takes input from makeAceObject
	#to calculate coverage for each base
	#output is a list
	
	tmp.pad.df <- data.frame(pad.pos=1:nchar(AceObject$padded.seq),
		base=strsplit(AceObject$padded.seq,"")[[1]],
		coverage=0,stringsAsFactors=FALSE)
	
	tmp.AF <- AceObject$AF[-grep('contig',AceObject$AF$read),]
	tmp.QA <- AceObject$QA[-grep('contig',AceObject$QA$read),]
	#tmp.read.len <- tmp.QA$stop.pos.clear + 1 - tmp.QA$start.pos.clear
	tmp.cont.stop <- tmp.AF$start.pos + tmp.QA$stop.pos.clear -1
	tmp.cont.start <- tmp.AF$start.pos + tmp.QA$start.pos.clear -1
	tmp.cont.start[tmp.cont.start < 1] <- 1
	
	if(length(tmp.cont.start)==0) return(NULL) #ie if no reads
	
	tmp.filter <- tmp.cont.stop > 1
	tmp.cont.stop <- tmp.cont.stop[tmp.filter]
	tmp.cont.start <- tmp.cont.start[tmp.filter]

	if(!identical(tmp.AF$read,tmp.QA$read)) error('AF does not match RD')
	
	for(i in 1:length(tmp.cont.stop)){
		#add 1 to coverage for each base for each overlapping read
		
		tmp.pad.df$coverage[tmp.cont.start[i]:tmp.cont.stop[i]] <- tmp.pad.df$coverage[tmp.cont.start[i]:tmp.cont.stop[i]] +1
		}
	list(coverage=tmp.pad.df[tmp.pad.df$base!="*",],padded.coverage=tmp.pad.df)
	}
#correlates with depth given by grouse1hc, but not equal to
#base positions seem fine, though if - on ref is out
#go for some heuristic approach..? 
#... eg only work on snps found where all pops coverage >= 10?
#Newbler must use the qual file as well...
#... doubt I have time to use that, just fudge it

